function f = LL_MTO2_BS(x,data)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%THIS FILE IS THE LIKELIHOOD FOR BOOTSTRAPPING THE ORDERED PROBIT MODEL IN 
%"Evidence of Neighborhood Effects from Moving to Opportunity: 
%     LATEs of Neighborhood Quality"
%by Dionissi Aliprantis and Francisca G.-C. Richter
%SEE THE INCLUDED readme.txt FILE FOR AN EXPLANATION OF THIS FILE 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%
%LABEL VARIABLES
%%%%%%%%%%%%%%%%
D10         =data(:,1);     %D \in {1, 2, ..., 10}
HHnoteens   =data(:,2);     %X_1 
HHnofam     =data(:,3);     %X_2 
HHvictim    =data(:,4);     %X_3 
basequal    =data(:,5);     %X_4 
sitedummy2  =data(:,6);     %BOSTON  SITE DUMMY
sitedummy3  =data(:,7);     %CHICAGO SITE DUMMY
sitedummy4  =data(:,8);     %LA      SITE DUMMY
sitedummy5  =data(:,9);     %NYC     SITE DUMMY
exp         =data(:,10);    %INDICATOR FOR MTO EXPERIMENTAL VOUCHER 
S8          =data(:,11);    %INDICATOR FOR SECTION 8 VOUCHER
weight      =data(:,12);    %WEIGHTS
mover       =data(:,13);    %INDICATOR FOR MOVE WITH VOUCHER

N=length(D10);              %NUMBER OF OBSERVATIONS
T=2;                        %NUMBER OF TYPES (t=1 for S8, t=2 for Experimental)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CREATING LATENT INDEX COMPONENT BY COMPONENT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%
%mu(X)
%%%%%%
bHHnoteens   =x(1);
bHHnofam     =x(2);
bHHvictim    =x(3);
bbasequal    =x(4);
bsitedummy2  =x(5);
bsitedummy3  =x(6);
bsitedummy4  =x(7);
bsitedummy5  =x(8);

mu=zeros(N,1);
mu(:)=bHHnoteens*HHnoteens(:) +bHHnofam*HHnofam(:) +bHHvictim*HHvictim(:) +bbasequal*basequal(:) ...
    +bsitedummy2*sitedummy2(:) +bsitedummy3*sitedummy3(:) +bsitedummy4*sitedummy4(:) +bsitedummy5*sitedummy5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COST FUNCTION (ie, CUTOFFS)
%
% !!!!!! OJO !!!!!!!!!!!!!
%
%HERE IT IS COST OF GOING FROM j-1 TO j
%IN PAPER IT IS COST OF GOING FROM j TO j+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bcutpoint(1)   =-1000000;
bcutpoint(2)   =x(9);
bcutpoint(3)   =x(10);
bcutpoint(4)   =x(11);
bcutpoint(5)   =x(12);
bcutpoint(6)   =x(13);
bcutpoint(7)   =x(14);
bcutpoint(8)   =x(15);
bcutpoint(9)   =x(16);
bcutpoint(10)  =x(17);

C=zeros(N,11);               %COST CUTOFFS
%should be index=-Inf --> C(:,1)=-Inf;
for j=1:10
    C(:,j)=bcutpoint(j);
end
C(:,11)=1000000;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COST-REDUCTION FROM SECTION 8 VOUCHERS (FOR HIGH-TYPES)
%
% !!!!!! OJO !!!!!!!!!!!!!
%
%HERE IT IS COST-REDUCTION OF GOING FROM j-1 TO j
%IN PAPER IT IS COST-REDUCTION  OF GOING FROM j TO j+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%g8cutpoint(1)    =x(18); - IRRELEVENT WITH NEW j INDEXING
g8cutpoint(2)    =x(18);
g8cutpoint(3)    =x(19);
g8cutpoint(4)    =x(20);
g8cutpoint(5)    =x(21);
g8cutpoint(6)    =x(22);
g8cutpoint(7)    =x(23);
g8cutpoint(8)    =x(24);
g8cutpoint(9)    =x(25);
g8cutpoint(10)   =x(26);

g8sitedummy2     =x(27);
g8sitedummy3     =x(28);
g8sitedummy4     =x(29);
g8sitedummy5     =x(30);

gamma8=zeros(N,11);          %EFFECT OF VOUCHER INCLUDING UNOBSERVED HETEROGENEITY
for j=2:10
    gamma8(:,j)=g8cutpoint(j) +g8sitedummy2.*sitedummy2(:) +g8sitedummy3.*sitedummy3(:) ...
                                           +g8sitedummy4.*sitedummy4(:) +g8sitedummy5.*sitedummy5(:);
end
gamma8(:,11)=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COST-REDUCTION FROM EXPERIMENTAL VOUCHERS (FOR HIGH-TYPES)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%g1cutpoint(1)   =x(31); - IRRELEVENT WITH NEW j INDEXING
g1cutpoint(2)   =x(31);
g1cutpoint(3)   =x(32);
g1cutpoint(4)   =x(33);
g1cutpoint(5)   =x(34);
g1cutpoint(6)   =x(35);
g1cutpoint(7)   =x(36);
g1cutpoint(8)   =x(37);
g1cutpoint(9)   =x(38);
g1cutpoint(10)  =x(39);

g1sitedummy2  =x(40);
g1sitedummy3  =x(41);
g1sitedummy4  =x(42);
g1sitedummy5  =x(43);

gamma1=zeros(N,11);          %EFFECT OF VOUCHER INCLUDING UNOBSERVED HETEROGENEITY
for j=2:10
    gamma1(:,j)=g1cutpoint(j) +g1sitedummy2.*sitedummy2(:) +g1sitedummy3.*sitedummy3(:) ...
                                           +g1sitedummy4.*sitedummy4(:) +g1sitedummy5.*sitedummy5(:);
end
gamma1(:,11)=0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MARGINAL BENEFIT (ie, LATENT INDEX) OF INDIVIDUAL n FROM MOVING FROM LEVEL j-1 TO LEVEL j 
%     (mu + gamma8_j^H*Z8 + gamma1_j^H*Z1 - C_j - V)
%
% !!!!!!!!! OJO !!!!!!!!!!!!!!!!!!!
%
%IN PAPER IT IS
%MARGINAL BENEFIT (ie, LATENT INDEX) OF INDIVIDUAL n FROM MOVING FROM LEVEL j TO LEVEL j+1 
%     (mu + gamma8_j^H*Z8 + gamma1_j^H*Z1 - C_j - V)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
index=zeros(N,11);          %LATENT INDEX FOR NBD QUALITY SELECTION
index(:,1)=-C(:,1);
for j=2:11
    index(:,j)=mu(:)+exp(:).*mover(:).*gamma1(:,j)+S8(:).*mover(:).*gamma8(:,j)-C(:,j);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CREATING TYPE PROBABILITIES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%
%SECTION 8
%%%%%%%%%%
b8HHnoteens   =x(44);
b8HHnofam     =x(45);
b8HHvictim    =x(46);
b8basequal    =x(47);
b8sitedummy2  =x(48);
b8sitedummy3  =x(49);
b8sitedummy4  =x(50);
b8sitedummy5  =x(51);

b8constant    =x(52);

index8=zeros(N,1);          %LATENT INDEX FOR SECTION 8 TYPE PROBABILITY
index8(:)=b8constant+b8HHnoteens.*HHnoteens(:) +b8HHnofam.*HHnofam(:) +b8HHvictim.*HHvictim(:) +b8basequal.*basequal(:) ...
           +b8sitedummy2.*sitedummy2(:) +b8sitedummy3.*sitedummy3(:) +b8sitedummy4.*sitedummy4(:) +b8sitedummy5.*sitedummy5(:);



%%%%%%%%%%%%%
%EXPERIMENTAL
%%%%%%%%%%%%%
b1HHnoteens   =x(53);
b1HHnofam     =x(54);
b1HHvictim    =x(55);
b1basequal    =x(56);
b1sitedummy2  =x(57);
b1sitedummy3  =x(58);
b1sitedummy4  =x(59);
b1sitedummy5  =x(60);

b1constant    =x(61);

index1=zeros(N,1);          %LATENT INDEX FOR EXPERIMENTAL TYPE PROBABILITY
index1(:)=b1constant+b1HHnoteens.*HHnoteens(:) +b1HHnofam.*HHnofam(:) +b1HHvictim.*HHvictim(:) +b1basequal.*basequal(:) ...
           +b1sitedummy2.*sitedummy2(:) +b1sitedummy3.*sitedummy3(:) +b1sitedummy4.*sitedummy4(:) +b1sitedummy5.*sitedummy5(:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CREATING LOG-LIKELIHOOD FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%
%COVARIANCE MATRIX "SIGMA"
%%%%%%%%%%%%%%%%%%%%%%%%%%
sigmaDD=1;
sigma88=1;
sigma11=1;

sigmaD8=x(62);
sigmaD1=x(63);

SIGMA8=zeros(2,2);
SIGMA8=[sigmaDD sigmaD8 ; 
        sigmaD8 sigma88  ; ];

SIGMA1=zeros(2,2);
SIGMA1=[sigmaDD sigmaD1 ; 
        sigmaD1 sigma11  ; ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PROBABILITY OF OBSERVATION n CONDITIONAL ON PARAMETERS FOR ALL j
%
% !!!!!!!!!!!  OJO  !!!!!!!!!!!!!!!!!!!!!!
%
%HERE
%Pr(D=j) = F_V(mu + gamma8_{j}^H*Z8 + gamma1_{j}^H*Z1 - C_{j}) 
%         - F_V(mu + gamma8_{j+1}^H*Z8 + gamma1_{j+1}^H*Z1 - C_{j+1}) 
%
%IN PAPER IT IS
%Pr(D=j) = F_V(mu + gamma8_{j-1}^H*Z8 + gamma1_{j-1}^H*Z1 - C_{j-1}) 
%         - F_V(mu + gamma8_{j}^H*Z8 + gamma1_{j}^H*Z1 - C_{j}) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%INFINITE BENEFIT FROM MOVING FROM D=0 TO D=1

Phij=zeros(N,1);
Phijplus=zeros(N,1);

    Phij(:)    =    Phij(:)+100000.*(D10(:)==1);
    Phijplus(:)=Phijplus(:)+index(:,2).*(D10(:)==1);
for j=2:9
    Phij(:)    =    Phij(:)+index(:,j).*(D10(:)==j);
    Phijplus(:)=Phijplus(:)+index(:,j+1).*(D10(:)==j);
end
    Phij(:)    =    Phij(:)+index(:,10).*(D10(:)==10);
    Phijplus(:)=Phijplus(:)-100000.*(D10(:)==10);

Pr=ones(N,1);
	Pr(:)=(normcdf(Phij(:)) ...
                      -normcdf(Phijplus(:))).*(S8(:)==0).*(exp(:)==0) ...
           + (mvncdf([Phij(:),  index8(:)], zeros(1,2), SIGMA8) ...
               -mvncdf([Phijplus(:),  index8(:)], zeros(1,2), SIGMA8)).*(S8(:)==1).*(mover(:)==1) ...
           + (mvncdf([Phij(:),  -index8(:)], zeros(1,2), SIGMA8) ...
               -mvncdf([Phijplus(:),  -index8(:)], zeros(1,2), SIGMA8)).*(S8(:)==1).*(mover(:)==0) ...
           + (mvncdf([Phij(:),  index1(:)], zeros(1,2), SIGMA1) ...
               -mvncdf([Phijplus(:),  index1(:)], zeros(1,2), SIGMA1)).*(exp(:)==1).*(mover(:)==1) ...
           + (mvncdf([Phij(:),  -index1(:)], zeros(1,2), SIGMA1) ...
               -mvncdf([Phijplus(:),  -index1(:)], zeros(1,2), SIGMA1)).*(exp(:)==1).*(mover(:)==0);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CONTRIBUTION TO LIKELIHOOD OF OBSERVATION n CONDITIONAL ON PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LL=zeros(N,1);              %LOG-LIKELIHOOD CONTRIBUTION OF OBSERVATION i
LL(:)=double(log(max(0.0000001,Pr(:)))).*weight(:);
             
f=double(-sum(LL(:)));



